function [fwhm, rpv, rpvmap] = smoothest3d(restc, msk)
% smoothest3d  - estimate smoothness map of residual time course data
%
% FORMAT:       [fwhm, fwhmmap, rpv, rpvmap] = smoothest3d(restc [, msk])
%
% Input fields:
%
%       restc      residuals time course data
%       msk        3D mask volume (logical)
%
% Output fields:
%
%       fwhm       global FWHM estimate
%       rpv        global resels per voxel estimate
%       rpvmap     resels per voxel map
%
% Note: code inspired by SPM2/5, see spm_est_smoothness.m

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   (~isa(restc, 'single') && ...
    ~isa(restc, 'double')) || ...
    ndims(restc) ~= 4
    error( ...
        'BVQXtools:BadArgument', ...
        'Function must be called with a 4-D single or double argument.' ...
    );
end
if nargin < 2 || ...
   ~islogical(msk) || ...
    ndims(msk) < 3 || ...
    size(msk, 1) ~= size(restc, 1) || ...
    size(msk, 2) ~= size(restc, 2) || ...
    size(msk, 3) ~= size(restc, 3)
    msk = [];
end

% get sizes
ts = size(restc);
vs = ts(1:3);
nt = ts(4);
nv = prod(vs);

% find mask
if isempty(msk)
    msk = 1:nv;
else
    msk = find(msk(:)');
end
nm = numel(msk);

% prepare summation arrays
vx = zeros(1, nm);
vy = zeros(1, nm);
vz = zeros(1, nm);
ss = zeros(1, nm);

% iterate over time points
for tc = 1:nt
    [dx, dy, dz] = gradient(restc(:, :, :, tc));
    dx = 2 * dx(:)';
    dy = 2 * dy(:)';
    dz = 2 * dz(:)';
    if nm == nv
        vx = vx + dx .* dx;
        vy = vy + dy .* dy;
        vz = vz + dz .* dz;
        ss = ss + restc(1+(tc-1)*nv:tc*nv) .^ 2;
    else
        vx = vx + dx(msk) .* dx(msk);
        vy = vy + dy(msk) .* dy(msk);
        vz = vz + dz(msk) .* dz(msk);
        ss = ss + restc(msk + (tc - 1) * nv) .^ 2;
    end
end

% normalize derivatives
vx = vx ./ ss;
vy = vy ./ ss;
vz = vz ./ ss;

% eliminate zero variance voxels
ina = isnan(vx) | isnan(vy) | isnan(vz);
vx(ina) = [];
vy(ina) = [];
vz(ina) = [];
msk(ina) = [];

% resels per voxel (resel) 
% resels = resels/voxel = 1/prod(FWHM)
% FWHM   = sqrt(4.ln2/|dv/dx|))
% fwhm   = 1/FWHM
%-----------------------------------------------------------------------
fwhm = sqrt([vx; vy; vz] ./ (4 * log(2)));
resel = prod(fwhm, 1);
rpvmap = zeros(vs);
rpvmap(msk) = resel;
iin = find(resel > 1);
rpvmap(msk(iin)) = 1;
fwhm(:, iin) = [];
resel(iin) = [];

% global equivalent FWHM {prod(1/FWHM) = (unbiased) RESEL estimator}
%-----------------------------------------------------------------------
fwhm = mean(fwhm, 2);
rpv = mean(resel);
fwhm = 1 ./ (fwhm * ((rpv / prod(fwhm)) .^ (1 / 3)));
